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1. INTRODUCTION 



Wildland fire is a complicated multiscale process. Fortunately, a practically important range of wildland fire behavior 
can be captured by the coupling of a mesoscale weather model with a simple 2D fire spread model ^Clark et alT] 
[1996a b). Weather has a major influence on wildfire behavior; in particular, wind plays a dominant role in the fire 
spread. Conversely, the fire influences the weather through the heat and vapor fluxes from burning hydrocarbons 
and evaporation. The buoyancy created by the heat from the fire can cause tornadic strength winds, and the wind 
and the moisture from the fire affect the atmosphere also away from the fire. It is well known that a large fire "creates 
its own weather." The correct wildland fire shape and progress result from the two-way interaction between the fire 
and the atmosphere ^Clark et al.,1996a|bt|2004(|Coen|2005/ . 

1.1. Origins and the current state of ttie modei 



WRF-Fire | |Mandel et a l."2009"i combines the Weather Research and Forecasting Model (WRF) fSk amarock et al. 



1008\ with a semi-empirical fire spread model, based on the level-set method. WRF-Fire has grown out of 
The NCAR's CAWFE code | |Clark et air]|1996 a b 2064| |Coen"]|2005| , which consists of the Clark-Hall mesoscale 



atmospheric model coupled with a tracer-based fire spread model, using the spread rate computed from McArthur's 
| |Noble et al.|1980i and later (Rothe rmel|1972 ) formula. Although the Clark-Hall model has many good properties, it 
is a legacy serial code, not supported, and difficult to modify or use with real data, while WRF is a parallel supported 
community code routinely used for re al runs. See (Coen and Patton 20101 for a further discussion of their relative 
merits in the wildland fire application. i jPatton and Coen,2004, i proposed a combination of WRF with the tracer-based 
model from CAWFE, formulated a road map, and made the important observation that the innermost domain of the 
weather code, which interacts directly with the fire model, needs to run in the Large Eddy Simulation (LES) mode. 
Patton ported the tracer-based code to Fortran 90, rewrote the heat flux insertion for WRF variables, and produced 
a prototype code coupled with WRF, which then served as the foundation of further development. However, instead 
of using the existing tracer-based code, the fire module in WRF-Fire was developed based on the level-set method 
fOsher and Fedkiw,|2003^ , among other reasons, because the level-set function can be manipulated more easily 
than tracers for data assimilation. Thus, only the fuel variables and the subroutine for the calculation of the fire 
spread rate remained from CAWFE. 

Since the paper | |Mandel et al.||2009) was written in 2007, a number of algorithms and other aspects have 



changed, new features were a dded, and WRF-Fire has been released as a part of WRF starting with version 3.2 
in April 2010 (Dudhia 2010; Wang et al. 20lOl. The latest version of WRF-Fire with new features and fixes which 
have not made it into the WRF download yet, plus additional visualization tools, guides, and diagnostic utilities, are 
available from the developers atlo penwfm. org| 

1 .2. Computationai simuiations 

While WRF-Fire takes advantage of the CAWFE experience, WRF is quite different from the Clark-Hall atmospheric 
model and the fireline propagation algorithm is also different. Thus, it needs to be demonstrated that WRF-Fire can 
deliver similar results as CAWFE, and WRF-Fire needs to be validated against real fires. | |Kim||2011 1 has shown 
that the level set method in the fire mod ule can advect the fire shape correctly, just like the CAWFE tracer code in 
dClark et al.|2004^ . i jJenkins et al.|2010) studied the role of wind profile on fire pr opagation speed and the shape of 
the fireline and demonstrated fireline fingering behavior on ideal examples, as in i jClark e t al. 1996a b). (Kochanski 
et al. 2010) compared simulatio n results with measurements on the FireFlux grass fire experiment (Clements et al." 



2007 ■ ( [Dobrinkova et al.|2010| simulated afire in Bulgarian mountains using real meteorological and geographical 
data, and ideal fuel data. (Beezley et al. 2010) simulated a fire in Colorado mountains using real data from online 
sources. A mesoscale simulation can run faster than real time on a small cluster | |Jordanov et al.|201 T ^ . 
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1.3. Related work 



Wildland fire models range from tools base d based on fir e spread rate formulas | |Rother mel"1 972' '1 983), such as 
BehavePlus ( |Andrews,2007) and FARSITE ( |Finney|1998) , suitable for operational forecasting, to sophisticated 3D 
computational fluid dynamics and combustion simulations suitable for research and reanalysis, such as FIRETEC 
| |Linn et a l. 20021 and WFDS (Mel l et al.| |2007). BehavePlus, the PC-based successor of the calculator-based 
BEHAVE, determines the fire spread rate at a single point from fuel and environmental data; FARSITE uses the fire 
spread rate to provide a 2D simulation on a PC; while FIRETEC and WFDS require a parallel supercomputer and 
run much slower than real time. 

The level set-method was used for a surface fire spread model in (iMallet et al. '20091 fFilippi et al.'2009i coupled 
the atmospheric model, Meso-nh, with fire propagation by tracers. Tiger (Mazzoleni and Giannino 2010) uses a 2D 
combustion model based on reaction-convection-diffusion equations and a convection model to emulate the effect 
of the fire on the wind. FIRESTAR | |Morvan and Dupuyl |2004) is a physically accurate wildland fire model in two 
dimensions, one horizontal and one vertical. UU LES-Fire (Sun et al. 2009) couples the University of Utah's Large 
Eddy Simulator with the tracer-based code from CAWFE. See the survey by | |Sullivan||2009^ for a number of other 
models. 



2. PHYSICAL FIRE MODEL AND FUELS 

The physical model consists of subroutine s computing the spread rate and the burn rat e, and it is essentially the 
same as a subset of CAWFE ^Clark et al.|[T996a,b)|2004;,Coen,2005) , | |Rothermel||1972^ , and BEHAVE ( Andrews] 
2007) . 



2.1. Fuel properties 

Fuel is characterized by a vector of quantities, which are given at every p oint of the domain. To simplify the 
specification of fuel properties, fuels are given as one of 13 ( |Anderson|1982| categories, which are preset vectors 
of values of the fuel properties above. These values are specified in an input text file, and they can be modified by 
the user. The user can also specify completely new, custom fuels. 



2.2. Fire spread rate 

Mathematically the fire model is posed in the horizontal {x,y) plane. The semi-empirical approach to fire 
propagation used here assumes that the fire spread rate is given by the modified Rothermel's formula 

5 = i?o (1 + + 0s) • (1) 
The spread rate depends on fuel properties, the wind speed U, and the terrain slope tan^, exactly as in fRothermel] 



1972>, except that some of the input quantities are metric so they are first converted from metric to English units 



(Ib-ft-min) to avoid changing the numerous constants in the computation from | |Rothermel 1972 ); the wind is limited 
to between to 30m/s; the slope is limited to nonnegative values; the fuel mass with moisture is given rather than 
dry fuel as in | |Rothermel|1972) ; and a new category is introduced for chaparral from | |Coen et al.|2001| eq. (1)). In 
either case, the spread rate can be written as 

S = max 1 5*0 , Rq + cmin {e, max {0,U}}'^ + dmax {0, tan (/)}| , (2) 

where So,Ro,b,c,d,e are fuel-dependent coefficients, which is how the spread rate is represented in WRF-Fire 
internally. These coefficients are stored for every grid point. At a fixed point on the fireline, denote by n the outside 
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normal to the fire region, if tine wind vector, and z the terrain height. The normal component of the wind vector, 
■ n, and the normal component of the terrain gradient, tan = Vz • n, are used to determine the spread rate, 
which is interpreted as the spread rate in the normal direction n. 

2.3. Fuel burned and heat released 

Each location starts with fuel fraction F = I. Once the fuel is ignited at a time U, the fuel fraction decreases 
exponentially, 

F(t)=cxpf-^^-^V t>U (3) 



where t is the time, ti is the ignition time, Fq is the initial amount of fuel, and T/ is the time constant of fuel, i.e., the 
number of seconds for the fuel to burn down to 1/e w 0.3689 of the original quantity. Since the fuel burns down to 
0.6 of the original quantity in 600s when w = 1000, we have 

(*-'.) igoo / {t-U) 

0.6 600 ™ = exp ' 



. Tf 
which gives 

_ 600 _ w 
^ ~ ~1000 In 0.6 ^ 0.8514' 

The input coefficient w is used in WRF-Fire rather than Tf for compatibility with existing fuel models and literature 
The average sensible heat flux density released in time interval (<, t + At) is computed as 

F{t)-F{t + M) 1 r. 2, , ... 

4>h^ — — wih (J/m /s) (4) 

^ At l + Mf ^ ' ^ ' 

and the average latent heat (i.e., moisture) flux density is given by 

F{t)-F{t + At) Mf+ 0.56^ ,,.2,, ,r. 

(h„ ^ — — T-^ — Lwi (J/m /s) (5) 

At 1 + A// 

where 0.56 is the estimated mass ratio of the water output from combustion to the dry fuel, and L = 2.5 • 10'' J/kg 
is the specific latent heat of condensation of water at °C, used for nominal conversion of moisture to heat. This 
computation is from CAWFE. 

3. MATHEMATICAL CORE OF THE FIRE MODEL 

The model maintains level set function ^, time of ignition tj, and the fuel fraction F. These quantities are all 
represented by their values at the centers of the fire mesh cells. 

3.1 . Fire propagation by the level set method 

This section follows i jMandel et al.|2009 |. Denote point on the surface by x = {x, y). The burning region at time t is 
represented by a level set function V = (x, t) as the set of all points x such that where V' (x, t) < 0. There is no 
fire at X if (x, t) > 0. The fireline is the set of all points x such that (x, t) = 0. Since on the fireline, the tangential 
component of the gradient V^p is zero, the outside normal vector at the fireline is 
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Now consider a point x {t) that moves with the fireline. Then the fire spread rate 5 at x in the direction of the normal 
n is 

^ ^ ax 



From tp (x (t) ,t) = and chain rule, 



d dip dip dx dtp dy 



dt' 



dip 



IIWI! n 



9x 



dip 



swvipw 



so, from and l|8), the level set function is governed by the partial differential equation 



dip_ 
~di 



+ 5(x)||V7M| =0, 



(7) 



(8) 



(9) 



called the level set equation | |Osher and Fedkiw|2003 |. The spread rate S is evaluated from ^ everywhere on the 
domain. Since S >Q, the level set function does not increase with time, and the fire area cannot decrease, which 
also helps with numerical stability by eliminating oscillations of ip in time. 

The level set equation is discretized on a rectangular grid with spacing (Aa;, Ay), called the fire grid. The level 
set function ip and the ignition time ti are represented by their values at the centers of the fire grid cells. This is 
consistent with the fuel data given in the center of each cell also. 

To advance the fire region in time, we use Heun's method (Runge-Kutta method of order 2), 



(10) 



^n+l/2 ^ _^ ^^p (-^n) 

= + At (^-F ir) + \f [r+"^)^ , 

The right-hand side F is a discretization of the term -5 \\Vip\\ with upwinding and artificial viscosity specifically, 

F{iP) = -5(1^ •7t,Vz • it) IjVV'll +eAiP, 

where it = VV'/I!VV'|| is computed by central differences and Vip = \y^ip,Vyip] is the upwinded finite difference 
approximation of Vip by Godunov's method |Osher and Fedkiw|2003| p. 58), 

W^i/j if ^^ip < and ^^ip < 0, 



V if V > and V^ip > 0, 
if v;^V < and V^?/' > 0, 



otherwise ip if 



(11) 



v>if v,^ < v:v 

where V+V and V~ip are the right and left one-sided differences 

ip{x + Ax,y) -ip{x,y) 



Ax 

Ip {x,y) -ip{x- Ax,y) 
Ax 
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and similarly for V+V and V'l/'- Further, 



is the gradient by central differences, e is the scale-free artificial viscosity (e = 0.4 in the computations here), and 

is the scaled five-point Laplacian of ^p. 

A numerically stable scheme that includes upwinding, such as Godunov's method l[Tl}, is required to compute 
IIVV^II in the level set equationj9}. However, it seems better to use standard central differences for VV' in the 
computation of the normal n in |[6]C which is needed to evaluate the normal component of the wind and the slope in 
<§■ 

Before computing the one-sided differences up to the boundary, the level set function is extrapolated to one layer 
of nodes beyond the boundary. However, the extrapolation is not allowed to decrease the value of the level set 
function under the value at either of the points extrapolated from. For example, when is the last node in the 
domain in the direction x, the extrapolation 

is used, and similarly in the other cases. This modification of the finite difference method serves to avoid numerical 
instabilities at the boundary. The extrapolation at the boundary effectively implements a free boundary condition. 
Without the stabilization, a decrease of tp at a boundary node, which often happens for nonhomogeneous fuels in 
real data, is amplified by extrapolation and because of upwinding, tp keeps decreasing at that boundary node forever 
and developing a large negative spike. 

The model does not support fire crossing the boundary of the domain. When t/' < is detected near the 
boundary, the simulation terminates. This is not a limitation in practice, because the fire should be well inside the 
domain anyway for a proper response of the atmosphere. 

The ignition time <j in the strip that the fire has moved over in one timestep is computed by linear interpolation 
from the level set function. Suppose that the point x is not burning at time t but is burning at time t + At, that is, 
'0(x, > and V'(x, t + At) < 0. The ignition time at x satisfies i){x,t,{x)) = 0. Approximating V' by a linear 
function in time, we have 

^(x,t,) - V>(x,t) _ ^P{:>^,t + At)-i;{K,ti) 
t, (x) - t ^ t + At- t, (x) 

and we take 

/ X 'ib(:x.,t)At 
tJx)=t+— , ^ ; — . 12 

3.2. Computation of the fuel fraction 

The fuel fraction is approximated over each fire mesh cell C by integrating i|3f over the fire region. Hence, the fuel 
fraction remaining in cell C at time t is given by 

1 ff l_,,Jj.^llM)dx. (13) 



area(C) J J Tf{x) 
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Once the fuel fraction is known, the heat fluxes are computed from Q and (5). This scheme has the advantage 
that the total heat released in the atmosphere after the fuel has completely burned is accurate, regardless of 
approximations in the computation of the integral |T3). We are looking for a scheme that is second order accurate 
when the whole cell is on fire, exact when no part of the cell C is on fire (namely, returning the value one), and 
provides a natural transition between these two cases. Just like the standard numerical schemes are exact for 
polynomials of a certain degree, the guiding principle here is that the scheme should be exact in a collection of 
(nontrivial) special cases. 

While the fuel time Tf can be interpolated as constant over the whole cell, the level set function V and the 
ignition time U must be interpolated more accurately to allow submesh representation of the burning area and 
gradual release of the heat as the fireline moves over the cell. Then, to compute the integral in |[T3), the cell C is 
split into 4 subcells C-, , and 

i;{x,t)<0 i/j(x,t)<0 

The level-set function tp is interpolated bilinearly to the vertices of the subcells Cj, and Tf is constant on each Cj. 
When the whole cell C is on fire (that is, V' < on all four vertices of C), ti is interpolated also linearly to the vertices 
of the subcells Cj. However, the case when the fireline crosses the cell C requires a special treatment of the ignition 
time til ti (x) has meaningful value only when the node x is on nodes on fire, -0 (x) < 0, and ti (x) on the fireline, 
i.e., when (x) = 0. Approximating both ^ and U in the fire region by a linear function suggests interpolating from 
the relation 

U-t = aP, (15) 

for some c. We interpolate on the grid lines between two nodes first. If both nodes are on fire, we interpolate 
ti bilinearly as before. However, when one cell center is on fire and one not, say ipi^i) > 0, ip{a2) < 0,. 



we find the proportionality constant c in ([15| from ti {a2) = cip{a2), and set <i (b) = cip{h) at the midpoint 
b = (ai + /2. In the case of interpolation to the node c = (ai + a2 + as + a.4) /4 between nodes ai, a2, aa, a4, 
we find the proportionality constant c by solving the least squares problem 



V(a,)<0 

and set again t^ (c) = cip (c). 

To compute the integral over a subcell Cj, we first estimate the fraction of the subcell that is burning, by 

area{(a;,y) e Cj : -0 (x,?/,t) < 0} _1( Efc=iV'(x/c)\ 



(16) 



where x^. are the the corners of the subcell Cj. The approximation is exact when no part of the subcell Q, is on 
fire, that is, all V (xfe) > and at least one V' (xfc) > 0; the whole Cj is on fire, that is, all V (xfe) < and at least one 
1/) (xfc) < 0; or the values tp (x^) define a linear function and the fireline crosses the subcell diagonally or it is aligned 
with one of the coordinate directions. 

Next, replace t (xk) by t when ^ (xfe) > (i.e., the node Xfe is not on fire), and compute the approximate fraction 
of the fuel burned as 
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This computation is a second-order quadrature formula when the whole cell is burning; it is exact when no 
part of the cell is burning; and it provides a natural transition between the two. Also, the calculation is accurate 
asymptotically when the fuel burns slowly and the approximation /? of the burning area is exact. 

Optionally, the fuel fraction remaining is computed first by approximating V and U and by linear functions using 
finite differences and then integrating analytically, first in the direction orthogonal to the fireline and then parallel to 
the fireline, thus splitting the fire region in the subcell into trapezoids. This method is exact when V and ti are linear 
functions, but it is more expensive. Accurate calculation of the fuel fraction and thus of the heat flux will be important 
for the simulation of crown fire, which is ignited when the ground heat flux exceeds a certain threshold value ^Clark] 
etal.|1996b} ). 



3.3. Ignition 

Typically, a fire starts much smaller than the fire mesh cell size, and both point and line ignition need to be supported. 
The previous ignition mechanism (Mandel et al. 2009) ignited everything within a given distance from the ignition 
line at once. This distance was required to be at least one or two mesh steps, so that the initial fire is visible on the 
fire mesh, and the fire propagation algorithm from Sec. |3.1| can catch on. This caused an unrealistically large initial 
heat flux and an accelerated ignition. 

The current ignition scheme achieves submesh resolution and zero-size ignition. A small initial fire is 
superimposed on the regular propagation mechanism, which then takes over. Drip-torch ignition is implemented 
as a collection of short ignition segments that grow at one end every time step. Multiple ignitions are supported. 

The model is initialized with no fire by choosing the level set function 7/;(x, to) = const > 0. Consider an initial 
fire that starts at time tg on a segment ab and propagates in all directions with an initial spread rate Sg until distance 
rg is reached. At the beginning of every time step t such that tg <t < tg + rg/Sg, we construct the level-set function 
of the initial fire, 

(x, t) ^ dist (x, a^)~Sg{t-tg) (18) 
and replace the level-set function of the model by 

ijj {x,t) mm{Tp (x,t) ,■0., (x,i)} . (19) 

For a drip-torch ignition starting from point a at time tg at velocity v until time th, the ignition line at time t is the 



segment a, a + v (min {t, th} - tg), and 18 becomes 



ipi (x, t) = dist [x, a, a + V (min {t, th} -tg)j - min {vg, Sg {t - tg)} 
followed again by |T9), at the beginning of every time step t such that 

tg <t <th+ Tg/Sg. 

The ignition time on newly ignited nodes is set to the arrival time of the fire at the spread rate Sg from the nearest 
point on the ignition segment. 

4. ATMOSPHERIC MODEL 

We summarize some background information about WRF from | |Skamarock et al.|[2008) to the extend needed to 
understand the coupling with the fire module. 
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4.1 . Variables and equations 



The model is formulated in terms of the hydrostatic pressure vertical coordinate 77, scaled and shifted so that ry = 1 
at the Earth surface and ry = at the top of the domain. The governing equations are a system of partial differential 
equations of the form 

(20) 

where $ = {U,V,W,(t)' ,Q, ^' ,Qrn)- The fundamental WRF variables are fi = n{x,y), the hydrostatic component 
of the pressure differential of dry air between the surface and the top of the domain, written in perturbation form 
= ji + fi', where ^ is a reference value in hydrostatic balance; U = ^u, where u = u{x,y,ri) is the Cartesian 
component of the wind velocity in the x-direction, and similarly V and W; O = ^i9, where 9 = 9{x,y,rj) is the 
potential temperature; cp = (t){x,y,ri) = <f) + ({>' is the geopotential; and Q,n = mm is the moisture contents of the 
air. The variables in the state $ evolved by l|20| are called prognostic variables. Other variables computed from 
them, such as the hydrostatic pressure p, the thermodynamic temperature T, and the height z, are called diagnostic 
varia bles. The variables that contain ^ are called coupled. The value of the right-hand side R ($) is called tendency. 
See (Skamarock et al. 12008 pp. 7-13) for details and the form of R. 

The system (|20| is discretized in time by the explicit 3rd order Runge-Kutta method 

$1 = + — 
3 

$2 = + 

$*+^* = $* + Ati? ($2) (21) 

where only the third Runge-Kutta step includes tendencies from physics packages, such as the fire module 
I jSkamarock et al.||2008[ p. 16). In order to avoid small time steps, the tendency in the third Runge-Kutta step 



also includes the effect of substeps to integrate acoustic modes. 



4.2. Surface sctiemes 

In real cases, non zero sf_sf ciay.physics should be selected to enable the surface model, allowing for 
proper interaction between the atmosphere and the land surface. In idealized cases, users have an option of 
the basic surface initialization, intended to be used without the surface model, or the full surface initialization 
(sf c_fuii_init=i). The latter allows for using all standard land surface models even in idealized cases. For 
idealized cases with full surface initialization, land surface properties like roughness length, albedo etc., are defined 
through the land use category The surface scheme utilizes a gridded array containing the number of landuse 
category, defined in a text file (landuse . tbl), which specifies the roughness length and other surface properties, 
both for the real and idealized cases. The land use categories may be also defined directly trough the namelist 
variables or read in from an external file containing a 2D landuse matrix (see also Sec.|7). 

5. COUPLING OF THE FIRE AND THE ATMOSPHERIC MODELS 

The terrain gradient is computed from the terrain height at the best available resolution and interpolated to the fire 
mesh in preprocessing. Interpolating the height and then computing the gradient would cause jumps in the gradient, 
which affect fire propagation, unless high-order interpolation is used. 

In each time step of the atmospheric model, the fire module is called from the third step of the Runge-Kutta 
method. First the wind is interpolated to a fixed height z/ above the terrain (currently, 6.1m following BEHAVE), 
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assuming the wind log profile 

const log — , z > zq, 



u{z) 



< z < zo, 



where z is the height above the terrain and zq is the roughness height. For a fixed horizonal location, denote by 
zi, ^2, ■ • ■ the heights of the centers of the atmospheric mesh cells; these are computed from the geopotential (jj, 
which is a part of the solution. The horizontal wind component U{zf ) is then found by interpolating the values 
U {zo) = 0,U (zi) , U {z2) , ... from the nodes logzo, logzi,logz2, ... to logz/; if zj < zo, we set U {zf) = 0. The 
y-component of the wind is interpolated in the same way. 
The fire model then makes one time step: 

1. If there are any active ignitions, the level-set function is updated and the ignition times of any newly ignited 
nodes are set following Sec. |3.3[ 



2. The numerical scheme l[To}-([lT| for the level set equation (9) is advanced to the next time step. 

3. The time of ignition set for any any nodes that were ignited during the time step, from ([12). 



4. The fuel fraction is updated following Sec. 3.2 



5. The sensible and latent heat flux densities are computed from Q and (5} in each fire model cell. 

The resulting heat flux densities are averaged over the fire cells that make up one atmosphere model cell and 
inserted into the atmospheric model, which then completes its own time step. The heat fluxes from the fire are 
inserted into the atmospheric model as forcing terms in the differential equations of the atmospheric model into a 
layer above the surface, with assumed exponential decay with altitude. Such a scheme is needed because WRF 
does not support flux boundary conditions. The sensible heat flux is inserted as the tendency of the potential 
temperature 0, equal to the vertical divergence of the heat flux, 

[x, y,z) = Re ($) H r — — exp 



dt ' ' ag{x,y,z) dz \ z^xt 



where Rq ($) is the component of the tendency in the atmospheric model equations (|20), a is the specific heat of 
the air, g{x,y,z) is the density, and Zext is the heat extinction depth, a parameter. The latent heat flux is inserted 
similarly into the tendency of the vapor concentration by 



d{mm) , ^ p , ^J,{x,y)(|)g{x,y) d 

[x, y, z) = Rq^ ($) H — — ^ — ^ — — exp 



z 



dt ' ' " Lg{x,y,z) dz \ z^xt 

where L is the specific latent heat of the air. 



6. SOFTWARE STRUCTURE 
6.1 . Grids and parallel execution 

Parallel computing imposes a significant constraint on user programming technique. WRF uses the RSL parallel 
infrastructure (Michalakes 20001. RSL divides the domain horizontally into patches. Each patch executes in a 
separate MPI process and is further divided into tiles, which execute in separate OpenMP threads. Communication 
between tiles is accomplished at the end of OpenMP parallel loop over tiles. The fire grid has refined tiles in the 
same location as atmospheric grid tiles. The patches are declared in memory with larger bounds than the patch size, 
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WRF: Runge-Kuttastep3 



Driver: get grid variables, get 
flags, interpolation calls, 
OpenlVIP loops, DM halos 



IVIodel: onetime step, winds 
in, lieat fluxes out 



Core: Advance the level set 
equation, ignition, compute 
fuel fraction. Dimensionless. 



Atmospheric pliysics: get 

temperatureand moisture 
tendencies from heat fluxes 



Fire physics: sensible and 
latent heat fluxes from fuel 
fraction, fire rate of spread 



Utilities: interpolation, WRF stubs, debug I/O,... 
WRF: error messages, log messages, constants,... 



Figure 1 : Structure of WRF-Fire software. All physics routines are in the dashed box. The Utilities layer is used by 
all other components above it. 



and communication between patches is accomplished by HALO calls (actually, includes of generated code), which 
update a layer of nodes beyond the patch boundary from other patches. The fire module computational code itself 
is designed to be tile-callable; it executes on a single tile, assuming it can safely read data from a layer of nodes 
beyond the tile boundary. The communication (OpenMP loops or HALO calls) occurs outside the computational 
routines; this means that whenever communication is necessary, the fire module must exit, and then continue from 
the correct code location on the next call. 



6.2. Software layers 

The fire module software is organized in several isolated layers (Fig.[l}. The driver layer serves to interpolate and 
otherwise translate the variables between WRF and the fire module, and it contains all exchange of data between 
the tiles in parallel execution. The rest of the code executes on a single tile, assuming that the needed values from 
neighboring tiles are already p resent. This structure is needed so that the rest of the code can conform to the 
WRF coding conventions (WRF Working Group 2|2007) . Only the driver layer depends on WRF; the rest of the fire 



module can be used as standalone code, independent of WRF WRF infrastructure is accessed only through stubs 
in the utility layer so that it can be easily emulated in the standalone code. The model layer is the entry point to the 
fire module. The core layer is the engine of the fire model, described in Sec. [3] The fire physics layer evaluates the 
fire spread rate and heat fluxes from fuel properties, and the atmospheric physics layer mediates the insertion of the 
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fire fluxes into the atmosptiere, as described in Sec. (5) One of ttie goals of the design is that the only components 
that need to be modified when the fire module is connected to another atmospheric model in future are the driver 
layer, the atmospheric physics layer, and WRF stubs in the utility layer. 



7. DATA INPUT 

WRF ideal run is used for simulations on artificial data. An additional executable, ideal . exe, is run first to create 
the WRF input. A different ideal . exe is built for each ideal case, and the user is expected to modify the source of 
such ideal case to run custom experiments. The ideal run for fire supports optional input of gridded arrays for land 
properties, such as terrain height and roughness height. This allows a user to execute simulations that go beyond 
what would normally be considered an ideal run and simplifies custom data input; the simulation of the FireFlux 
experiment was done in this way |Kochanski et al.|20i0| . 

In a real run, the WRF Preprocessing System (WPS) | |Wang et ar]|2010 ! Chapter 3) takes meteorological and 
land-use data in a number of commonly used formats and prepares it for WRF to use as initial and boundary 
conditions. WPS has been extended to process fine-scale land data for use with the fire model, such as topography 
and fuel. 
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